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ABSTRACT 

We describe a technique for solving for the orbital elements of multiple planets from radial velocity 
(RV) and/or astrometric data taken with 1 m/s and /xas precision, appropriate for efforts to detect 
Earth-massed planets in their stars' habitable zones, such as NASA's proposed Space Interferometry 
Mission. We include details of calculating analytic derivatives for use in the Levenberg-Marquardt 
(LM) algorithm for the problems of fitting RV and astrometric data separately and jointly. 

We also explicate the general method of separating the linear and nonlinear components of a model 
fit in the context of an LM fit, show how explicit derivatives can be calculated in such a model, and 
demonstrate the speed up and convergence improvements of such a scheme in the case of a five-planet 
fit to published radial velocity data for 55 Cnc. 

Subject headings: astrometry — planetary systems — methods: data analysis — methods: numerical 
— techniques: radial velocities 



1. INTRODUCTION 



1.1. 



Fitting Kepelerian Curves 

The disco very of over 27 mul tiple-planet systems in re- 
cent years (jWright et al.l l2009h has required algorithms 
for disentangling the radial-velocity signature of such 
complex systems. Because the parameters describing a 
radial velocity (RV) or astrometric curve are nonlinear, 
there is no way to fit for them analytically, and they 
must be found through an algorithmic search. Fitting 
for the Keplerian parameters of a single orbital compan- 
ion is usually straightforward given a good period guess, 
and if necessary a "brute force" mapping of the \ 2 space 
is usually not computationally prohibitive. Fitting mul- 
tiple planets involves searching a correspondingly higher- 
dimensional space and can require substantial computing 
time. 

There is an art to searching such x 2 spaces efficiently, 
and in this context there are many "tricks" for finding 
the global minimum. F or instance, a Lomb-Scargle peri- 
odogram (|Scarglejll982[ ) is often used to identify promis- 
ing periods for prospective planets, and all of the tallest 
peaks can be used as starting guesses for the fitting algo- 
rithm. In hierarchical systems, the dominant planet can 
be fit for alone, its signal subtracted from the data, and 
additional planets can be searched for among the resid- 
uals. This process can then be repeated until all of the 
planets have been identified, and then a full, multi-planet 
fit on the original data sta rting at the values found for 
the individual planets, fe.g. iFischer et al.ll2008l ). 

In many multi-planet systems, planet-planet interac- 
tions can significantly alter the radial velocity (RV) or 
astrometric signature of the system. In such cases where 
interactions are important, a full dynamical (Newtonian) 

1 Townes Fellow, Space Sciences Laboratory, UC Berkeley 



fit involving an n-body code must be used to properly fit 
the data and to ensure the short- and long-term stabil- 
ity of the solution. Even in these cases, a multi-planet 
Keplerian (kinematic) fit, which simply adds the reflex 
signatures of single planets and ignores planet-planet in- 
teractions, is still useful for efficiently identifying planets 
and providing good initial guesses to the n-b qdy codes. 

The Levenberg-Mar quardt method (LM; iLevenberel 
fl94l lMarauardtl ll963) is an efficient algorithm for find- 
ing a local minimum in a nonlinear x 2 space (given good 
guesses), and is well-sui ted for the app lication of RV 
and astrometric fitting (|Press eTaI][l99^ . As a con- 
crete example, we refer in this work to a useful IDL 2 
impleme ntation of this tec hnique, MPFIT by Craig Mark- 
wardt 3 (|Markwardtl 12003 ) . Like most implementations 
of the algorithm, MPFIT requires a user-defined function 
that accepts, as an argument, trial values for the pa- 
rameters being solved for, evaluates the model at those 
values, and returns the corresponding residuals to the 
data. MPFIT then uses this information to step through 
parameter space and locate the minimum in x 2 using a 
combination of Newton's method and a steepest-descent 
method. Uncertainties in these parameters can then be 
calcu lated by mapping t he x 2 space near this minimum 
(e.g. I Wright et al.l l2007f) or through error "bootstrap- 
ping " (e.g. lButlere t al. 200j). 

The user-defined function in MPFIT also optionally 
returns the values of the derivatives of these residu- 
als with respect to the parameters being fit, computed 
"analytically" or "explicitly" . Absent these derivatives, 
MPFIT will calculate numerical derivatives with small 

2 IDL is a commercial programming language and environment 
by ITT Visual Information Solutions, http://www.ittvis.com/idl/ 

3 Available at http://purl.com/net/mpfit. MPFIT is a port 
of MINIPACK-1 from FORTRAN, and is also available in C and 
Python. 
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steps in the fitted parameters and calculating the result- 
ing change in the residuals, he user must choose step sizes 
that are not too large — thereby missing fine structure 
in the x 2 space — or too small — increasing compute 
time and potentially losing numerical precision. Adding 
explicit-derivative capability to MPFIT, or any implemen- 
tation of the LM algorithm, obviates the need for explic- 
itly providing step-sizes, and can greatly improve per- 
formance in terms of the number of steps taken and the 
CPU time consumed per step. 

The LM algorithm is useful for finding the best-fitting 
parameters of a nonlinear model to a set of data from 
a "frequentist" perspective. A Bayesean approach can 
provide for more robust estimates of parameter uncer- 
tainties, especially when those uncertainties are large. 
One Bayesean method of exploring complex or high- 
dimensional spac es is the Markov chain Monte Carl o 
method fMCMC: [Metropolis et al.1H95llHastingslll970D . 
which has been pro ductively employed in the context of 
orbital fitting (e.g Ifb7dl [200l iDriscoll fe Fischerl l200l 
iBalan fc LahavH2008D . 

Future space missions, such as NASA's Space Interfer- 
ometry Mission (SIM Lite) , will obtain /xas astrometry of 
nearby stars, referenced to an inertial astrometric grid. 
These measurements will be sufficiently precise to detect 
Earth-mass planets with orbital periods shorter than the 
mission lifetime. Since multiple-planet systems are com- 
mon, 4 interpreting these data in conjunction with precise 
radial velocity data will require algorithms that can effi- 
ciently and robustly search the large nonlinear parameter 
space of multiple-planet systems. 

In this paper, we describe the method of efficiently 
fitting multi-Keplerian models to such high precision 
RV and astrometric data separating the parameters 
into linear and nonlinear components in the context 
of a LM algorithm, and provide the explicit deriva- 
tives used in such a fit. With modifications, the prin- 
ciples here_caii_jls^_bj_j£pjied to M CMC methods, as 
well (jCatanzarite. Zhai. fc Shaoll2009t iBakos et al.ll2009L 
e.g.). 

1.2. Plan 

We begin with an elementary exposition to familiarize 
the reader with our notation (detailed in Table [2|) and 
provide context for the later discussion. In §SJ we expli- 
cate the method of exploiting linear parameters in the 
Kepler problem using the example of an RV time series. 
The calculation of explicit derivatives in W2A\ is general 
to any application of the LM algorithm where the model 
contains both linear and nonlinear parameters. We ap- 
ply this method to the problem of astrometric data in £|3] 
and specifically to the problem of combined astrometric 
and RV data in £13.31 We discuss nonlinear terms relevant 
to /ias astrometric work and how to accommodate them 
in §3] In f}5]we quantify the improvement in speed and 
convergence from exploiting linear parameters and from 
the use of explicit derivatives in the LM algorithm. 

2. AN EXAMPLE: RADIAL VELOCITIES 

2.1. Costs and Benefits of Exploiting Linear 
Parameters 

4 At least 1/4 of known planetary sy stems show evidence of 
multiple companions (Wright ct al. 2007) 



When fitting RV data, there are 5n + 1 Keplerian pa- 
rameters to be fit, where n is the number of planets: 5 P, 
the period of the planet's orbit; K, the semi-amplitude of 
the radial velocity signal; e, the eccentricity of the orbit; 
lo, the argument of periastron; t p , a date of periastron 
passage; and 7, the apparent radial velocity of the center 
of mass of the system 6 

Exploiting linear parameters, as described below, re- 
duces the search space to 3n dimensions (corresponding 
to P, e, and t p for each planet) and combines the other 
orbital elements into a set of linear parameters, which 
can be solved analytically (and therefore quickly and ex- 
actly) at each step in the search. The exact, analytic 
solution of these linear parameters greatly increases the 
speed and stability of search algorithms, but at a cost: 
the nonlinear parameters cannot be varied independently 
of the linear parameters. 

Further reduction in the number of nonlinear pa- 
rameters per planet is certainly possible. By exploit- 
ing an epicyclic or harmonic series expansion one can 
reduce the problem to only one n o nlinea r pa rame- 
ter pe r planet, Pj. ICumming et alj (|2003l ) and iFordl 
(2008) analyze the RV problem in the case of a cir- 
cular orbits, and discuss the relationship between pe- 
riodograms and Bayesian ap proaches to orbit fitting. 
iKonacki fe Macieiewskil (|1999f ) pursued a method analo- 
gous to ours in their approach to RV curv e fitting, and 
IKonacki. Macieiewski. fc Wolszczanl (|2002D did the same 
for astrometry. Such approximations offer a different set 
of costs and benefits to the one presented here. For in- 
stance, when K and e are both large, it may require a 
large number of terms to adequately describe a set of 
RV data. We may pursue such an approach in a future 
version of our code. 

The linear parameters in our treatment are not coeffi- 
cients in a series expansion; rather, we recast the prob- 
lem, separating linear parameters, which can be solved 
for exactly with linear algebra, from nonlinear param- 
eters, which must be solved for algorithmically (using 
LM). As described in §2.2, the linear parameters are al- 
gebraic combinations of K, ui, 7, and an optional trend 
parameter, while P, tp, and e arc nonlinear parameters. 

Two complications are introduced by exploiting lin- 
ear parameters. First, because the linear parameters are 
computed analytically, and not algorithmically, their er- 
rors and covariances with the nonlinear parameters are 
not computed automatically in the procedure outlined 
here. (They also cannot be held fixed in the fit, although 
the issue of fixing the trend parameter can be finessed - 
see §2.6). The second complication is that computation 
of the explicit derivatives for the LM algorithm is not 
straightforward: in the context of the algorithm, the Ke- 
plerian model is not a simple function of the 3n nonlinear 
parameters being fit, but also depends on the data. For 
example, increasing the nonlinear parameter e slightly 
not only changes the model because it is more eccentric, 
but also because at this new value of e the linear pa- 
rameters have different values. These complications can 

5 Occasionally, an additional "trend" parameter is used to fit out 
long-term RV trends caused by massive, long-period companions. 
The other two elements, i, the inclination, and f2, the position angle 
of the ascending node, can only be determined astrometrically. 

e In practice, 7 is degenerate with an arbitrary instrument- 
dependent RV offset. 
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be avoided if one uses the linear basis of the orbits pre- 
sented here in an ordinary, nonlinear 5n + 1 parameter 
fit, but the convergence and speed-up benefits will not 
be as great. 

The sections below describe a general method of calcu- 
lating an n-planet Keplerian radial velocity model given 
values for the 3n nonlinear parameters, and of calculating 
the derivatives of that model with respect to those pa- 
rameters. These equations can be used in a user-defined 
function for input to an LM minimization code, or a 
brute-force x 2 map. 

2.2. Linear Parameters in the Radial Velocity Problem 
We wish to find the parameters of the n-planet model 

n 

u(t) = [Kj (cos^j + fj (£)) + ej coswj)] +j + d- (t-to) 

(1) 

which best fits the set of observed radial velocities v, mea- 
sured at times t, with uncertainties a in a least-squares 
sense. Here, Kj and uij are the usual Keplerian parame- 
ters for planet j, 7 is the time- independent velocity off- 
set, d is the trend parameter, fj(t) is the true anomaly of 
planet j at time t, and to is a conveniently chosen epoch 
of the observations. The true anomaly is defined im- 
plicitly in terms of the other three Keplerian parameters 
{Pj,t v j, and ej) through the relations 



tan 



' 1 + e j 
1 — e-t 



■ tan 



Ej(t) 



(2) 



Ej(t) -ej sin Ej(t) 



2ir(t - t Pi . 
Pi 



Mj(t) (3) 



In Eq. [3] (known as Kepler's Equation) Ej is called the 
eccentric anomaly of planet j, and M is known as the 
mean anomaly. 7 
We identify the linear parameters by rewriting Eq. [1] 

as 

n 

u(t) = ^2[hjcosf j (t) + cjsmfj(t)}+v + d-(t-t ) (4) 

j'=i 



where 



and 



hj = Kj cos £ 



Cj = —Kj sinwj 



vq = 7 + Kj ej cos uij 
3=1 



(5) 
(6) 

(7) 



These linear parameters can be converted back to Ke- 
plerian orbital elements through the relations 



Kj = ,/h*j + c* 



tanwi = 



(8) 
(9) 



7 Instead of t p , many authors (especially dynamicists) prefer 
to parameterize orbits in terms of the mean longitude at epoch, 
defined as M(to) + Q + w. 



(where ujj is chosen so that sinwj has the sign of the 
numerator) and 



7 = vo 



Kj6j COS Wj 



(10) 



The masses of the orbiting planets can be inferred from 
their corresponding semi-amplitudes Kj, defined for a 
single-planet system with planet mass m and stellar mass 
m* as 

r3 = 2nG ( m 3 sin 3 ? \ 
P(l- e 2)l \(m*+m) 2 ) 

where G is Newton's gravitational constant. The frac- 
tion in parentheses is known as the mass function of the 
system. 

The problem can now be divided into two parts: an 
algorithmic search through parameter space for the best- 
fit nonlinear parameters P. 



j ■ 



z.j , and t Pt j with a computer 



routine such as an LM or an MCMC code, and at each 
step in that search an analytic solution for the linear 
parameters that fit best there. 

2.3. Solving for the Linear Parameters 

Given a set of values for the nonlinear parameters, we 
can fit for the linear parameters in Eq. 0] through \ 2 min- 
imization. We denote the row vector of linear parameters 



(3 = {/ii,ci,/i 2 ,c 2 
We define \ 2 t ne usual way: 

JV 



■ h„,, c„, v , d} 



X" 



E 



(Vk - u{t k )y 



(12) 



(13) 



and minimize it with respect to each of the linear param 
eters in fj simultaneously: 

dx 2 _ 2 Vk ~ u{t k ) du 
dpi ~ f-1 o% dfji 



= 



(14) 



We can express this more compactly by invoking matrix 
algebra. For the problem of Keplerian orbits, we define 
the matrix F as 

cos/ M cos/i j2 .. 
sin/i,! sin/1,2 •■ 
cos/2,1 cos/2,2 ■■ 
sin/2,1 sin/2,2 •■ 



cos/„,i 
sin/„,i 
1 

t\ — to 



COS /„, 2 
sin f n .2 

1 

h — to 



(15) 



where fj k = fjit^)- This allows us to write the model 
velocities at times t (Eq. [J) as 

it=f3F (16) 

We also define the diagonal weight matrix W such that 

W k i = Ski/o-l (17) 

where we have used the Kronecker delta symbol. We can 
then write the system of equations in Eq. [T?] as 

dxfl 

dfj 



-2(?-/?F)WF T = d 



(18) 
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Solving for /3 we have 

(3 = «WF T e (19) 

where we have denoted the error matrix (also called the 
variance-covariance matrix) 

e = (FWF T )-' (20) 

Eq. Q1J] represents the general solution to the linear 
least-squares problem for an appropriately defined F. In 
the context of Keplerian fits, given a set of nonlinear 
parameters P, e, and t p for each planet, the remaining 
Keplerian orbital elements can be found from (3 using 
Eqs.EHEIl 

2.4. Calculation of Explicit (Analytic) Derivatives for 
Use in the LM Algorithm 

The derivative of the model velocities u with respect 
to any nonlinear parameter can be found from Eq. 1161 



F + /3 



dF 



(21) 



d.u 

dx dx dx 
where x stands for any of the nonlinear parameters (here, 



3 ■ 



, or £p 



From Eq. [19] we have 



dp 
dx 



vW 




(22) 



From the definition of a matrix inverse we can express 
the last derivative as 



de 

dx V dx 



(23) 



Using Eq. [20] we then have 

^ = -e4- (FWF T )e 

dx dx 



(24) 



Eqs. fT6H24l are not specific to the RV Kepler problem, but 
are a general method of calculating explicit derivatives 
in a model with both linear and nonlinear parameters, 
and so can be applied to any analogous problem. For 
instance, the problem of fitting an orbit using astrometric 
data also has linear parameters as we show in Ej3j 

2.5. Explicit Derivatives for the Radial Velocity Model 

For the case of a Keplerian RV model, from Eq. [T5l we 
have 



dF 

dx 



-sin/i.i/^j 

COS/l,l/( a 
-sin/2,1/2,1 
COS/ 2 ,l/2,l 

■ sin 







-sin/i i2 / li2 

COS /i )2 i 1 , 2 

-sin/ 2i2 /2,2 
003/2,2/2,2 

sin/ n , 2 /^ 2 

COS/„ j2 /4,2 






where /j k = dfj.k/dx. Note that since /j k 



(26) 



refers to the 



parameter of a different planet (e.g., dfj^/dPi — when 
j =/= I). This means that the matrix dF/dx has only 
two nonzero rows. We can therefore suppress subscripts 
below for clarity. 
We can calculate the nonzero derivatives as 



^ = dl + d£dE 
dx dx dE dx 

where, from Kepler's Equation (Eq. [3]), we have 



dE_ 
dP 



dE 
dt p 

dE 
~de~ 



-2n(t-t p )/P 2 
1 — e cos E 

-2n/P 
1 — e cos E 

smE 



1 — e cos E 



and from Eq. [2] we have 

df__ 
dE 



1 + e 1 + cos / 
1 - e 1 + cos E 

From Eq. [2] we can also write 

df _ 2tan(£/2)cos 2 (//2) 



de 



(l-e)vT 



(27) 

(28) 
(29) 

(30) 

(31) 
(32) 



but it is more computationally convenient to note that 
this happens to simplifiy to 



df df smE 
~de = ~d~E\-e 2 



Finally, 



df_ = df_ 
dP dU 



= 



(33) 



(34) 



true anomaly of planet j, it vanishes when x refers to a 



Eqs. |2"2TIM1 can be used to calculate the terms in Eq.l2~ll 
yielding the explicit derivatives used by LM method fit- 
ting routines, such as MPFIT. 

2.6. Variations on the n-Planet RV Model 

In the context of the Kepler problem, the above equa- 
tions include a parameter for a linear trend in the data. 
In practice, fitting for such a trend will only occasionally 
be necessary. When not needed the d parameter in (3 and 
the bottom rows of the F and dF j dx matrices can simply 
be left out of the calculations. Similarly, the trend pa- 
rameter can effectively be held fixed at a nonzero value 
by simply subtracting the desired value from the data 
before fitting. 

These matrices can also be easily extended to handle 
the case of combining data from multiple telescopes be- 
tween which there exist RV offsets. This is accomplished 
by extending the data vectors v, t, and a to include data 
from all telescopes, and extending the vector j3 to include 
a separate offset parameter for each telescope after the 
first. The corresponding rows of F must then be filled 
with l's in those columns corresponding to data from 
the appropriate telescope, and 0's elsewhere. Naturally, 
in dF/dx the elements of these rows are all 0. 
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3. APPLICATIONS TO ASTROMETRY 
3.1. Astrometry Alone 

The above method can also be extended to the prob- 
lem of fitting for the Keplerian elements of an orbit from 
astrometric data. Here we present a simplified model 
of /xas astrometric data of the sort that might be pro- 
vided by a space observatory such as NASA's Space In- 
terferometry Mission (SIM Lite). We anticipate mod- 
ifying our algorithms for a more realisti c model of in- 
terferometric data and its n o ise sources dSozzettil [20051: 
Eriksson fc Lindegrenl l2007t ICatanzarite. Law, fc Shad 
2008) in the near future. 

The linear basis for astro metric coordinate s are the 
Thiele-Innes constants (e.g. iBinne ndiu^ll96l . and arc 
well docu mented as useful tools for astrometric curve- 
fitting (e.g. ICasertano et al.ll2008h . The astrometric per- 
turbations caused by an orbiting companion can be de- 
scribed in terms of six astrometric orbital elements: in 
addition to e and t p , we have a, the semi-major axis of 
the star's apparent orbit on the sky in units of arc; f2, 
the longitude of the ascending (approaching) node (mea- 
sured as a position angle on the sky); i, the inclination 
of the orbit on the sky (such that i = corresponds 
to a face-on, clockwise orbit); and o>», the longitude of 
periastron of the star's orbit. 8 

The Thiele-Innes constants are defined in terms four 
of the astrometric elements of the star's orbit about the 
secondary: 



A = a{ cos COSW, — sin Q sin a;* cos i) (35) 

B — a{ sin f2 cos a>* -I- cos SI sinw, cos i) (36) 

F = a(— cos J! sin a;* — sin fl cosw* cos i) (37) 

G — a{~ sin f2 sin w„ + cos Q cos cos i) (38) 

C = a sin uj* sin i (39) 

H = a cos lu* sin i (40) 

These constants can be quickly computed using rotation 
matrices: 

A B C 

F G H 

asinzsinfi — asinicosfi acosz 



where R is the 3-D rotation matrix 
cos £7 sin ft 

R z {tt) 



= aR z {u*)R x {€)R z ({l) 
(41) 



sinf2 




cos il 




(42) 



and W* is the argument of periastron of the orbit of the 
star. 

We can transform the Thiele-Innes constants back to 
Keplerian orbital elements of the planet with the rela- 
tions: 

~Q Tp 

tan(w* + fi) = (43) 



tan(w» — Q) - 



A + G 
-{B + F) 



A-G 

( _ , _ {A - G) cos(^ + n) 
U " 2 J ~ {A + G) cos(w* — n) 



(44) 
(45) 



8 The orbital parameters of the star and those of the unseen 
companion are all identical except a, which differs by a factor of 
m/m„, and cj, which differs by 7r. 



and 



a= (A cos lu* — F sin lu*) cosfi— 
(A sin oj* + F cos ) sin Q sec i 



a;* 



(46) 



(47) 



where the quadrants of — f2 and +Q are determined 
by the signs of the numerators in Eqs. [43] & HU These 
relations leave a ±7r ambiguity inw,,w, and f2 that can 
only be resolved by radial velocities, without which con- 
vention dictates that we choose the solution with Q < it. 

The Thiele-Innes constants C and H are closely related 
to the c and h constants of Eqs. [5] & [5] The set of six 
constants have the identity 



A + B 



C 2 



F + G 



H 1 



(48) 



In astrometry there are five parameters that describe 
a star's motion in the absence of orbiting companions: 
Aaocosj and ASq, the difference between the true and 
nominal position of the system at to; fi a and (x$, the 
proper motions in the RA and Dec directions; and zu, the 
parallax of the system. Our model for the astrometric 
displacement of a star due to parallax, proper motion, 
and a system of unseen planets in terms of the Thiele- 
Innes constants at times f is: 



= E]=MjX jtk + F j Y j!k }+ 
A8 + zuILs,k + M<5( r fc - to) 



(49) 



Aa k cosS = J2"=ii B j x j,k + Gj Y jM+ (50) 
Aa cos S + zuU a ^ k + [i a (r k - to) 

where X and Y are the so-called elliptical rectangular 
coordinates, defined as 



X j. k = cos Ej(r k ) 



Y j,k = \ 1 ~ei sin Ej(T k ) 



(51) 
(52) 



where E is the eccentric anomaly and the quantities ILj,^ 
and II^^. refer to the astrometric displacements due to 
parallax in the a and 8 directions 9 during observation 
k, which are given by [cite Supplement to Astronomical 
Almanac here] : 



II a ,fe = r x (r k ) sin a - r y (T k ) cos< 



(53) 



ILjfc = (r x (rfe) cosa + r y (r k ) sina) sin 6 - r z (r k ) cos<5 

(54) 

Here (r x ,r y ,r z ) represent the Cartesian components in 
equatorial coordinates of the position of the observatory, 
r, at time r with respect to the Solar System barycenter 
(in units of AU when w is in arcsec). These values for 
the Earth are available from the NASA Jet Propulsion 
Laboratory Solar System ephimerides, 10 but for /zas and 
spaceborne work the precise position of the observatory 
itself is required. 

Note that since a is the apparent semi-major axis of the 
star's orbit in units of arc, its relationship to the mass 
of the secondary depends on the method of astrometry 
used. For astrometric perturbations due to an unseen 
planet, (that is, absolute astrometric displacements with 

9 In this paper the bare symbols a and S will always refer to the 
nominal right ascension and declination of a star at t$ , absent the 
effects of parallax and astrometric displacement from companions. 

10 http: //ssd. jpl.nasa.gov 
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respect to the sidereal frame, as measured by SIM Lite) 
we have from Kepler's Third Law 



(m* + m) ; 



(55) 



where a is measured in arcseconds when P is measured 
in years, w is the parallax in arcseconds, and m is the 
mass of the unseen companion and m* is the mass of the 
star in solar massses. 11 

We can now extend Eq[T|5] to the case of 2-D data by 
defining our vector of N measurements taken at times t: 

[A<5i, A<$2 ■ ■ ■ A<5at, Aai cos 5, Act2 cosS . . . Accn cos 6] 
and our model with linear parameters for n planets: 



(3 — [Ax, Bi, Fx, G\ . . . A n , B n , F n , G n 
ASq , Aa cos S, us, Ha, ro ] 



(57) 



and matrix F 
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3.2. Astrometry in Arbitrary Coordinates 

Astrometry does not always deliver contemporaneous 
(AacosJ, AS) pairs at a common time t. In the general 
case, a baseline determines the 1-D displacement of a star 
from some reference at an intermediate position angle on 
the sky (i.e., not necessarily [as in the case for AS] or 
7r/2 [as for Aacos^]). We can combine Eqs.09]&[5O]to a 
more general form to accommodate a heterogeneous set 
of such data : 

pe.k = YJj=i [(AjXj,k + FjYj,k) cos 9 k + 
{B J X^ k + G j Y ]k )fiind k } + 

(AS +wU s ,h + M«( T fc- *o)) cos 6 k+ 
(Aa cos S + wli ajk + ^a{T k - to)) sin6> fc 

(65) 

where p is the separation and 9 the position angle of the 
measurement. 

Interestingly, we could also use the other two Thiele- 
Innes constants to achieve the same result by defining 
a new linear parameter scheme where the astrometric 
displacements are described as: 

Pe,k = Ylj=i[HjSj,k + CjTj,k]+ 

(A5 + mlls,k + M(5( T fc - *o)) cos#k+ 
(Aa cos 5 + mll a ,k + Ma( T fc _ *o)) sin6» fe 

(66) 

where S and T are defined for the kth measurement and 
jth planet as 

T 3,k 



(58) 

The nonzero components of dF/dx can be calculated 
from: 



dX 

dP 
dX 

dt p 

dX 

de 
dY 

dP 
dY 

dt p 

dY 

de 



dE . 
: — — sin E 
dP 



dE_ 
dtp 
dE 
de 



sinE 



■ sin E — 1 



\/l — e 2 cos E 



= yl — e 2 cos E 



= yl — e 2 cos E 



dE 

dP 
dE 

dtp 

dE 



e sin E 



de vT 



(59) 
(60) 

(61) 
(62) 
(63) 
(64) 



and Eqs. 



11 In the case where the binary orbit is measured as a sep- 
aration and position angle of one star with respect to another 
(i.e., relative astrometry) the measured separation is given by 
a 3 = ro 3 (m* +m)P 2 . The application of the techniques here for 
multiple-planet systems with relative astrometry is not straightfor- 
ward, and is beyond the scope of this manuscript. 



sin(f^ 



- 9k) csci; 
9k) cot ij 



sin(Qj 
cos(f2j 



9k ) cot ij 
0k) esc ij 



Xj,k 

Yj,k 



(67) 

This scheme uses only two linear parameters, instead of 
four, per planet, and so will not be as efficient, but it 
is still useful because it can be easily combined with the 
radial velocity scheme of H2.2I We demonstrate this in 
^3.31 with a procedure that can accommodate any com- 
bination of radial velocity and astrometric data. 

We can recover the parameters a and ui from the linear 
parameters C and H given the nonlinear parameter i by: 



sin 2 ii 



tan uii 



-Cj_ 
-Hi 



(68) 
(69) 



where ojj is chosen so sinWj has the same sign as — Cj. 

3.3. Combining Astrometry with Radial Velocities 

Combining astrometric data with radial velocity 
data is not as simple as combining the RV-only and 
astrometry-only schemes outlined above, because the six 
linear parameters (A, B, F, G, and c and h) are a com- 
bination of only five Keplerian elements (K, a, f2, u>, and 
i), and the problem is thus overconstrained. One solu- 
tion would be to minimize \ 2 subject to the appropriate 
constraints using Lagrange multipliers, but the resulting 
set of nonlinear equations may not be guaranteed to have 
a unique solution and would be difficult to solve for the 
case of an arbitrary number of planets. 
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Another solution is to adapt the C and H constants to 
accommodate general astrometric data, and these con- 
stants are closely related to the scheme for RV data used 
in $2.2\ To do this, we define our vector of measurements 
to be 

[v,Pe\ (70) 

where the velocities are taken at times t and the astrom- 
etry at times f , and define our model, as in Eq. 1161 



(3 — [Hi, Ci, H-2, C 2 ■ ■ ■ H n , C n , vq, d, 
ASq , Aa cos 5, n$,ii a , zu] 

From Eqs. El M HB & E3 we have: 

hj = -XjHj 



(71) 



(72) 
(73) 



metric measurements, read 

Si,i 
Ti,i 
£2,1 
T 2 ,i 



T n .l 



S\ t 2 
Tl,2 

02,2 
T22 






cos 9i 
sin 0i 

(ri — to) cos bi (T2 — to ) cos c 2 

(ti - t ) sin 9i (t 2 - t ) sin 2 

ILj,i cos 61 + n Qj i sin 61 Us,2 cos 62 + n Qj2 sin( 



S n ,2 
T n ,2 




cos 82 
sin 82 
(t2 — io) cos 1 
(t 2 - t ) sin ( 



For the nonlinear parameters Pj , T p 
nonzero rows of dF/dx can be calculated from 



j, and ej, 



(76) 
the two 



where we have introduced a A, a combination of nonlinear 
orbital parameters of planet j which has units of velocity 
and has the value 



2vrAU 



(74) 



when the estimated parallax, w, is expressed in arcsec- 
onds. 

The parameters Hj, Cj, and vq can be transformed 



into mh j (m 



h and 7 with Eqs.[EHH&[I2H71 



The appearance of w in the definition of Xj (and there- 
fore in Cj and Hj) indicates a fundamental nonlinearity 
in the combined RV-astrometry problem: the parallax 
is not a truly linear parameter. However, because it is 
nearly linear, if a good estimate of the parallax is avail- 
able then this can be used in the formula for Xj (we indi- 
cate the approximate nature of the parallax term with a 
tilde). Once the parallax is solved for more precisely, the 
fit can be re-run with an improved estimate of w. This 
procedure should converge very quickly, and we will use 
it again to deal with other, smaller nonlinear terms in 

m 

The first columns of the F matrix, corresponding to 
the radial velocity measurements, now read 



-cos/i,iAi 
sm/i,iAi 

-COS/24A2 
sin/ 2 ,iA 2 



cos/ n ,iA r , 
sm/ n ,iA„ 
1 

ti — to 



Q 





- cos/i >2 Ai 
sin/i j2 Ai 

- cos/ 2i2 A 2 
sin/ 2i2 A 2 

- cos /„ )2 A n 
sin/ nj2 A„ 

1 

t2 — to 









(75) 



and the rest of the columns, corresponding to the astro- 



d_ 

dx 
d 

dx 
where 



: cos/(i)A) 
(-sin/(i)A) 



■sin/(t)/'(t)A + cos/(t)^ (77) 
ax 

■cos/(f)/'(t)A-sin/(i)^ (78) 
dx 



dX 

dP 
dX 

de 



A 
P 

eX 



1 



dX 
dtt 



dX _dX _ 
di dt-n 



(79) 
(80) 
(81) 



and 



dS/dx 
dT/dx 

cos(f2 
— sin(fi - 



) esc 1 
) cot i 



cos(r2 



! cot 1 

) esc i 



dX/dx 
dY/dx 



(82) 

and Eqs. |2~TH3"41 [B^HMl Now we have introduced two 
additional nonlinear parameters, f2 and i. 



ds/dn 
dT/dn 



— sin(fi — 8) esc i — cos(51 

— cos(f2 — 8) cot i sin(il 



) cot i 
) esc i 





' X ' 




Y 



(83) 



and 



dS/di 
dT/di 



cos(fl — 8) esc i cot i sin(fi — 8) esc 2 i 
sin(fi — 8) esc 2 i cos(J7 — 8) esc i cot i 



X 
Y 



(84) 



4. NONLINEAR TERMS 

4.1. Sources of Nonlinearity 

There are several small nonlinear terms which are im- 
portant at the m/s and /zas level, especially for the 
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nearby and high-proper motion stars likely to be ob- 
served by SIM Lite. 

Secular acceleration — A star with significant proper 
motion will have the radial component of its space ve- 
locity change with position on the sky, resulting in a 
secular change in the radial velocity up to ~1 m/s/yr for 
the most extreme cases. Secular acceleration is given, to 
first order, by: 

v r = D\i 2 (85) 

where, v r is the bulk radial velocity of the star, 12 D is 
the star's distance 13 and \i is the total proper motion 
in radians per unit time. This term will be absorbed 
into the linear parameter d, if present, and so could be 
ignored. 

Parallax changes — The change in parallax of nearby 
stars due to their radial velocity may be of order 0.3 
/uas/yr. The shape of the parallactic motion is also a 
function of position on the sky, and thus of the proper 
motion. These changes can be of order 5 /xas/yr. The 
radial velocity term is given by: 



D 



(86) 



Proper motion changes — The flip side of secular ac- 
celeration is proper motion change due to change in dis- 
tance. This effect can be of order 3 /ias/yr and is given 
by: 



l l = ~H 



D 



(87) 



Curvilinear effects — The curvilinear nature of spher- 
ical coordinates can produce what are essentially nonlin- 
ear terms depending on how astrometric displacements 
are defined. These effects are on the same order as the 
above proper motion changes, and more pronounced near 
the poles. For /zas astrometry, it suffices to handle these 
effects by employing a rectilinear grid: 

In this work, astrometric displacements labeled 
Aacos<5 and AS do not strictly refer to changes in the 
right ascension and declination of the star, but refer to 
displacements along rectilinear axes along those dimen- 
sions at the nominal position of the star. That is, if the 
unit vector pointing to the nominal position of the star 
from the Solar System barycenter is defined: 

p = [cos a cos 8, sin a cos 8, sin 8] (88) 

then the unit vectors pointing east and north are given 

by 

a=[0,0,l]xp (89) 
S=pxa (90) 
and the astrometric displacements are given by 



AS = p' -8 
Aa cos8 = p' ■ a 



(91) 
(92) 



12 This quantity v r defines the true radial velocity of the system 
with respect to the Solar System barycenter. It differs from the 
spectroscopic parameter 7 in that the latter is often measured with 
respect to a fiducial frame and can include non-Doppler effects 
such as instrumental offsets, gravitational redshift, and convective 
blueshift. 

13 D is distinguished here from the inverse parallax to -1 simply 
for convenience of units. 



where p' is the displaced position of the star and S and 
a are constant. 

Inter feromteric cross terms — The 1-D interferometric 
measurement of astrometric displacement on the sky may 
be complicated by the motion of the reference and target 
stars. That is, the calculation of 9 in Eq. 1651 may require 
proper motion advanced or parallax corrected positions 
in a manner specific to the details of a particular instru- 
ment's measurement of 0. These cross terms are likely 
to be small, and so they can estimated and refined in the 
same manner as the other nonlinear terms, if necessary. 

RV- astrometry cross terms — The RV semi-amplitude 
K is related to the astrometric semi-major axis, a, by 
the parallax, and so w is not strictly a linear parameter 
in the combined astrometry-RV problem. This effect can 
be large if the parallax is small or not known, and should 
not be ignored for any system. 

Relativistic Terms — Gravitational deflection by Solar 
System objects and relativistic stellar aberration produce 
large, time dependent astrometric displacements that de- 
pend on the position of the star, and are therefore slightly 
nonlinear. Because these displacements can be calcu- 
lated to better than /zas precision given an estimate of 
the star's true position to arcsecond precision, these ef- 
fects are ignored here. 

4.2. Incorporating Nonlinear Terms 

Implementing these small nonlinear effects in our 
model requires good estimates of the astrometric param- 
eters (when astrometric data is available, these parame- 
ters can be estimated from a first-pass solution assuming 
no planetary companions). The system is then solved us- 
ing these estimates to calculate the second-order terms 
above. Below, we indicate these estimated astrometric 
parameters with a tilde to distinguish them from the 
solved parameters. These estimates can then be iter- 
atively refined if necessary, but convergence should be 
very fast for SIM Lite data. 

We can include these nonlinear terms by making the 
following substitutions to the bottom 3 rows of the ma- 
trix F in Eqs. [55] & US 



T fc - to -»• 1 - (r fc - io)-k (n - t ) (93) 



2D 



n 



ot,k ' 



1 - (rk - ioW ) n ajfe 



l-(7*-*,)g 



n<5 h 



(94) 



(95) 



and in the definitions of H k (Eqs. [53] fc I5~4"]) : 

n a ,fc^II Q , fe + (Aa cosS + (r k - t )fi a )r-p (96) 

ILj, fe ->ILj,fc + ( ASq +(T k -t )fls)f-p (97) 
where the quantities in parentheses have units of radians. 
The estimated terms Aao cos 8 and ASo will likely be zero 
at first, but may be iteratively refined with the other 
estimated parameters. 

Finally, the secular acceleration can be accommodated 
by simply subtracting off the appropriate, approximate 
linear trend to the RV data, or, if the trend parameter is 
present, allowing it to be absorbed in d. 

This procedure of estimating and refining linear astro- 
metric terms will also work for the nonlinear term intro- 
duced by the appearance of w in H and C in Eqs. [72] & 
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[72l Alternatively, the parallax can be treated as a non- 
linear parameter from the outset and eliminated from the 
linear coefficients in (5 entirely. 

5. IMPLEMENTATION 

5.1. Public code 

We have implemented these algorithms into a set of 
IDL software routines which we have made available pub- 
licly for the fitting of radial velocity data. 14 The soft- 
ware package includes RVLIN, the user-defined function 
that may be passed to MPFIT, and RV_FIT_MP, a "wrap- 
per" routine that employs MPFIT to fit a multi-planet 
Keplerian model to a user-supplied set of radial veloc- 
ity data. We anticipate maintaining and improving this 
package, and eventually incorporating astrometric and 
transit data analysis. This code, or components of it, 
are currently used by members of the California Planet 
Search. Below, we discuss our software package's perfor- 
mance. 

5.2. Speed-up From the Use of Explicit Derivatives 

The use of explicit derivatives above in an LM code 
significantly speeds up the algorithm by avoiding unnec- 
essary calculation of numerical derivatives. To quantify 
this improvement, we tested two cases with published RV 
data, the 2-planet system HD 21 7107 dVotrt et al.l|2Q 05D 
and the 5-planet system 55 Cnc (Fi scher et al.ll2008f l (in 
both cases assuming no jitter). We combined the Lick 
and Keck RV data sets for both cases, solving for all 
5n + 2 parameters (including the RV offset between the 
two telescopes). We ran over 2,500 trials, where in each 
trial we started the search with different initial guesses 
for the orbital elements, each randomly drawn from a 
normal distribution with a width given by the uncer- 
tainty in each parameter and centered on its best-fit 
value. 15 We compared the total time taken for these 
trials on a 2.6 GHz Intel Core 2 Duo MacBook Pro run- 
ning IDL 7.0 using explicit derivatives to the time taken 
on the same machine with same initial guesses using nu- 
merical derivatives. 

In the case of HD 217107, the use of explicit derivatives 
sped up the calculation by a factor of 2.3. The improve- 
ment in the case of 55 Cnc was even larger, a factor of 
4. The total number of steps taken by MPFIT to converge 
on a solution was similar in the cases with and without 
explicit derivatives, indicating that our step sizes for the 
numeric derivatives were well chosen. 

5.3. Convergence Benefits of Exploiting Linear 
Parameters 

We employed a custom version of the multi-planet 
fitting routine often employed by the California and 
Carnegie Planet search as a baseline to test the improve- 
ment in convergence of the algorithm described in this 
w ork. This routine, whic h de rives from code described 
inlMarcv k. Butlerl (|1992l ) and lValenti. Butler, fc Marcvl 
(|1995f ). is essentially an LM algorithm for searching all 

14 Available at http://exoplanets.org/code/ 

15 Although dynamical fits are more precise, here we are only 
concerned with the algorithm's convergence in the region of the 
global x 2 minimum for a purely Keplerian fi t. T he initial guesses 
for the parameters in this and the test in § 15.31 were thus drawn 
near their values at this (presumed) global minimum. 



5n + 1 parameters with carefully chosen step-sizes for 
numerical derivatives. 16 

Usin g the same ha rdware described above, we again 
fit the IFischer et all (|2008l ) 55 Cnc RV data from two 
telescopes. We drew the initial guesses for every param- 
eter from normal distributions centered on the Keple- 
rian best-fit values and with width given by sa x , where 
a x is t he parameter uncertainty quoted in IFischer et all 
(2008) and s is a scale factor. We varied s smoothly 
from (where the initial guesses were the best-fit values 
exactly) to 10 (where every parameter is independently, 
randomly altered by lQa x ) over 4000 trials. 

We did not attempt to fit for the telescope offset at 
each trial, as this particular routine was not optimized 
for such a task. We also did not vary the 7 parameter 
from its best-fit value (though we did fit for i t), as no 
uncertainties were given for 7 in IFischer et ail (|2008l ). 

We compared the results of this custom code to our 
LM code employing linear parameters ( §5.1|> with both 
explicit and numerical derivatives. In some sense, this 
is not a fair test since the linear parameter models do 
not require or accept initial guesses for K and u>, and 
thus the initial guesses are closer to the best-fit values 
(since they differ in only 15 dimensions space, whereas 
the guesses for the nonlinear routine they differ in 25 di- 
mensions.) We thus also ran an additional test, where 
the fully-nonlinear routine was provided the best-fit val- 
ues for K and uj for each planet, and thus had only 15 
parameters varied (but still had all 26 parameters to fit). 
With this exception, we provided the same initial guesses 
to each of the four schemes in each trial. 

We recorded the final x 2 reported by each routine for 
each trial and compared this with Xmin> the best-fit value. 
We deemed any trial for which (x 2 — Xmin) < 2 to be a 
"successful" convergence on the correct parameters. 17 At 
each tested value of s, we weighted the full set of trials 
by a Gaussian window 75 trials wide (corresponding to 
~ 0.19c) and calculated the fraction on trials that were 
successful. 

In Figure Q] we plot the fraction of successful tests as 
a function of s. The use of numerical derivatives had no 
overall effect on the convergence of the linear pramater 
routine, but was slower, on average, by a factor of 4. 

Even with initial guesses IOct^ from the best-fit values 
in all 15 parameters, the linear parameter routine found 
the global minimum of \ 2 m roughly half of all trials. 
The fully-nonlinear routine, searching a 26-dimensional 
space, required much better initial guesses to achieve con- 
vergence. With the same 15 parameters varied, fits with 
initial parameters off by 3o~ x had only a 20% chance of 
properly converging. With all 25 orbital parameters var- 
ied, the initial guesses had to be within l.b-o~ x of their 
proper value to have a 50% chance of convergence. With 
26 dimensions to search, there are many wrong paths 
for the LM algorithm to follow away from the global \ 2 
minimum. 

6. CONCLUSIONS 

16 That is to say, the only significant diff erences between the 
baseline code and the code described in i) l5.1l is the exploitation of 
linear parameters. 

17 ^min ~ 2910. Our results are only weakly sensitive to the 
precise definition of "successful" 
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Fig. 1. — Sensitivity of RV-only fitting algorithms to ini- 
tial guesses on the 5-planet, two-telescope RV data for 55 Cnc 
lIFiscner et al. 1 120081) . Initial guesses to four fitting routines were 
randomly varied from their best-fit values by various factors of their 
respective uncertainties. The linear parameter fitting algorithm 
described in this work converged on the best solution for a wide 
range of initial guesses, including ~ 50% of cases with guesses 10<r 
from nominal. A full, 26-parameter nonlinear fit required guesses 
within 1.5-cr of the best-fit value, or within 2a if only the 15 truly 
nonlinear parameters were varied. 



a full, 26-parameter search requires all 26 parameters to 
be specified within 1.5-cr x . 

The principle improvement from use of linear parame- 
ters comes from reducing the search space for an n-planet 
model. In the case of RV fitting, the reduction is from 
exploiting the linearity of 2n + 1 of the 5n + 1 fitted pa- 
rameters, leaving only 3n nonlinear parameters to be fit 
algorithmically. The problem of fitting astrometric or- 
bits can be similarly treated by exploiting the linearity 
of 4n + 5 of the 7n + 5 parameters, leaving, again, only 
3n nonlinear parameters to be with a nonlinear fitting 
routine. When combining RV and astrometric data, as 
will be the case for SIM Lite, only 2n + 6 of the In + 6 
model parameters are usefully linear under the scheme 
described here, leaving 5n nonlinear parameters. 

We have derived a general expression for the explicit 
model derivatives for models employing linear parame- 
ters in nonlinear fits, appropriate for application in the 
LM method. This result is general and can be applied to 
problems beyond Keplerian fitting - indeed to any model 
with both linear and nonlinear parameters. 



Applying linear parameters to the problem of fit- 
ting Keplerian curves to radial velocity and astrometry 
data significantly improves the efficiency and reliability 
of multi-planet Keplerian fitting routines. This tech- 
nique can be applied to many methods of searching this 
complex, nonlinear, multi-parameter \ 2 space, including 
the Levenberg-Marquardt method, Markov chain Monte 
Carlo algorithms, and brute force approaches. Table [1] 
summarizes the various schemes used in this work. 

We have identified the nonlinear terms relevant for /zas 
astrometry and m/s radial velocity work (such as that 
by SIM Lite and its supporting RV data), and shown 
how to incorporate these terms into a linear parameter 
scheme. In £12.41 we have provided analytic forms for ex- 
plicit derivatives relevant to various applications of the 
Kepler problem. 

In the case of RV-only data, use of explicit derivatives 
can speed up a fitting routine by a factor of 2-4, de- 
pending on the number of planets being fit. Use of linear 
parameters greatly improves the convergence properties 
of a multi-planet fitting routine. In the case of an ac- 
tual 5-planet fit, a linear parameter model requires ini- 
tial guesses within 10a x of their correct values in only 15 
parameters to have a 50% chance of convergence, while 
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TABLE 1 
Summary of Linearization Schemes 







Radial Velocities 


Radial Velocities 


Astrometry 


&: Astrometry 



Inputs v(t),fi,,z&, Aa(f),A(S(f),pQ,(J5,i) r ,ro v(t), Pg{f), /Is, v r , to 

Nonlinear parameters Pj,t p j,ej Pj,t p j,ej Pj,t p j,ej,Qj,ij 



Transformed linear parameters Cj,hj, vo — » oJj,Kj,y Aj, Bj , Fj , Gj — ► , , , aj Cj , Hj , vq —> ojj , ( m + m .yi > 7 

Linear parameters d A<5p , Aojq , fig, fi a ,zu d, ASp , Aqp , fig , fi a , to 

Note. — The transformed linear parameters vq and 7 and all of the linear parameters appear once per system. The subscript 
j on the other parameters indicates that there are n such parameters, one for each companion in the system. 
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TABLE 2 

Variables and symbols used in this manuscript 



Symbol 


Meaning 


Example equation 




A tilde indicates a first approximation as opposed to a fitted parameter 




9.3 




A, B, C, F, G, H 


Thiele-Innes constants 




11 




a 


Astrometric semi-major axis of a star's orbit in units of arc 




55 




a 


Nominal right ascension of a system at the epoch of observations, to 








a 


Constant unit east vector 


ED 





Vector of linear parameters 


El 


c 


Linear parameter in RV fitting corresponding to a component of a planet's RV signature 








d 


Linear parameter in RV fitting corresponding to an RV trend 








D 


Distance to a system 


EH 


Act cos <5, A8 


Measured astrometric displacements in units of arc in the ct and 8 directions 


cm a 


so: 


.5 


Nominal declination of a system at the epoch of observations, to 








8 


Constant unit north vector 


M 


Ski 


Kronecker delta 


El 


7 


Constant (and often instrument dependent) offset in a set of RV data 




m 




E 


Eccentric anomaly of a planet (a function of time) 








e 


Eccentricity of a planet 




i 




F 


Matrix defined such that the model u = /3F 


EI& 


58 


f 


True anomaly of a planet (a function of time) 




m 




e 


Position angle of astrometric displacement such that 9 = refers to A<5 


12 


h 


Linear parameter in RV fitting corresponding to a component of a planet's RV signature 




m 




i 


Inclination of a planet's orbit with respect to the sky 


mi 


j 


Subscript indicating a quantity corresponds to the fth planet 




i 




K 


RV semi-amplitude of a star's orbit due to planet 




i 




k 


Subscript indicating a quantity corresponds to the kth observation 


E] 


A 


A combination of nonlinear orbital parameters 


EH 


M 


Mean anomaly 


1 


771* 


Mass of the primary component of a binary or planetary system 




55 




m 


Mass of a smaller component of a binary or planetary system 




55 






Proper motion. The total proper motion is given by = /i^ -1- ^ 




58 






Number of planets in a system 




1 




iV 


Number of observations being fit 


E] 


y 2 
X 


1 lll_ oL'dlJOLJL 


El 


p 


Period of a planet 




3 






TTnit" n nm i n :. 1 nn^ifinn \rr i r t tnv nf n Qf"nT" in rmv^rr'ftni'rir pfiiiatdnal cn-irriin 

l_l 111 L 11U111111C11 UUCl L1U11 V11UU1 Ul < I oldl 111 UCll yUUllllllj LU UCILU1 Idl LUU1 U11K1HS 




88 




lie. , 115 


Hi i n !■ ."i nil t: ( ax ."inn" 1 arin tm n^r^TiKitnj null" upriillpf^ir"' mdtifiri 

1 HllL 11U11-5 1 Ul llllll_ ClllU 1 I Ul_.5ljl 1 Ulllti 11111 1 UCll C111CIU 111, HlUtlUll 


53 


& 


HI 




The mathematical constant 










t-^avtillav nr n drtiiom 
1 dl (lllelX Ul ' 1 bySUclll 


.55 


& 


58. 


1 KE ,. . 1 L 2S V / 


The 3-D rotation matrix about the x- or is-axis. 




12 




' J * X) ' y ) ' Z 


(1 ritiov^/nt'OVA/ Tinni ."inn in nanri^fiTitTif' fun n ."n."!:. 1 ^norHiripi + oc; I n fnnr'ri.in Mr hinpl 

V 1 L1SL1 VdLUI V UUol L1U11 111 UCll yL.l_llllll_ 1UUHLUI 1C11 LUU1 U111CI LI,.- I CI 111111U1U11 Ml I lllll_ 1 


53 


& 
65 


nn 


P6 


A/I qii ttti nctrntnrtrir rlit;nlpr , mTir>Ti1" ill tno HinT"Tinn nf nncitinn uniylp H 

IVlldDUl 111 doll U1I1L Ul 11 Uli5LllClH_llll_.ll I 111 1111. 1111 l_l_ L1U11 Ul UUolllUll CllltlH. 1/ 








Scale of random deviation of initial guesses from nominal in units of Ux 




5.3 


q rp 


Astrometric terms in the F matrix 




67 






Ivlcasurcment uncertainty 


g 


& 


17 


(Jx 


Uncertainty in an orbital parameter, x 


u 


& 


17 


% 


Time of a radial velocity observation 




1 




til 


Fiducial time at the epoch of the observations 


93 


tr, 


Time of periastron passage of a planet 




m 




r 


Time of an astrometric observation 


\M\M 


u 


Model values, such as velocities or astrometric displacements 


m 


&E1 


1; 


Measured radial velocities 


EJ 


T'O 


Linear parameter corresponding to a constant offset in a set of RV data 




1 




V r 


Radial velocity of a system's barycenter with respect to Solar System barycenter 




85 




W 


Diagonal matrix containing weights of measured data 




17 




x,y 


Elliptical rectangular coordinates of a planet (functions of time) 


51 


& 


52 




As a variable, can stand for any parameter, such as Pj, f Pl j, or ej 




21 




n 


Position angle of the ascending (approaching) node of a planet 




11 




a; 


Argument of periastron of a planet's orbit 




1 






Argument of periastron of a star's orbit due to a planet, u) = u>* + tt 





